
function out = frameSegmentationPDN( str )
global HD HR
HD = frameFilterD('linear');
HR = frameFilterR('linear');
%get input image, renormalized to [0,1]
I = readImage( str );
f = I;
psf = fspecial('gaussian',2,2);
u0 = f;
[M N] = size(f);
f  = imfilter(f, psf);
f = (f-min(f(:)))/(max(f(:))-min(f(:)));
usmooth =f;

figure(1);
imshow(I,[]);
%set parameter values.
global nPhases nLevels
global M N
nPhases = 4;
nLevels = 1;     %level of framelet transform
nIterations = 50;
tao = 0.4;
sig = 0.4;
gamma = 4;
beta =0.5;
nu = 0.1;
alpha = 0.5;

%  c = kmeansC(f,nPhases);
% c= [0,1];
%  c=[   0    0.4478    0.6927    0.8460];
%  c = sort(Kmeans(f,4)); 
c=[0.1744    0.4086    0.6160    0.8250];
u = cell(1,nPhases);
d = cell(1,nPhases);
for i=1:nPhases
    u{i}=f;
    d{i}= frameDec2(usmooth,nLevels,HD);
end

%set gw, the edge indicator function.
tm1 = d{1}{1};tm1(1:M,1:N) = 0;
sigma = 1;
gw = zeros(M,N);

for k = 1 : nLevels
    tmp = d{1}{k};
    tmp(1:M,1:N)=0;
    for i = 1 : numel(HD)
        for j = 1 : numel(HD)
            gw = gw+abs(tmp((i-1)*M+1:i*M,(j-1)*N+1:j*N)).^beta;
        end
    end
end

gw = 1./(sigma*gw+1);
gw = (gw - min(gw(:)))./(max(gw(:))-min(gw(:)));
gw = repmat(gw,numel(HD),numel(HD));
ubar = u;
uold = u;

y = cell(1,nPhases);
for i = 1 : nPhases
    for j = 1 : nLevels
        y{i}{j}= zeros(M*numel(HD),N*numel(HD));
    end
end

wu = cell(1,nPhases);
Eold = 0;
Ecur = 0;

%BEGIN ITERATION
%-----------------------------------------------------------
for k = 1 : nIterations
    k
    tic
    for i = 1 : nPhases
        wu{i} = frameDec2(ubar{i},nLevels,HD);
        y{i} = cellPlus(y{i},cellMultiply(wu{i},sig.*gw));
    end
    Ecur = frameEnergy(u,wu,gw,c,f,nu);
    ETF(k) = Ecur;
    err(k)=abs(Ecur-Eold)/Eold;
    err(k)
%          if err(k) <1e-6
%              break;
%          end
    Eold = Ecur;
    for i = 1 : nPhases
        for j = 1 : nLevels
            s = zeros(M,N);
            tmpy = y{i}{j};
            for l = 1 : numel(HD)
                for q = 1 : numel(HD)
                    if (l ~=1 || q ~=1 )
                        s = s + tmpy((l-1)*M+1:l*M,(q-1)*N+1:q*N).^2;
                    end
                end
            end
            s = max(sqrt(s),1);
            s = repmat(s,[numel(HD),numel(HD)]);
            s(1:M,1:N) = max(tmpy(1:M,1:N),1);
            y{i}{j}=y{i}{j}./s;
        end
    end



W=[];
for i = 1 : nPhases
    ry = frameRec2(cellMultiply(y{i}, gw), nLevels, 0, HR);
    u{i} = u{i} -tao*nu*(c(i)-f).^2 - tao*ry;
    W = [W;u{i}(:)'];
end
Wp = frameProj(W);

for i = 1 : nPhases
    u{i} = reshape(Wp(i,:),[M,N]);
end

theta = 1/sqrt(1+2*gamma*tao);
tao = theta*tao;
sig = sig/theta;
ubar = u;
%     W = [];
%     for i = 1 : nPhases
%         W = [W;u{i}(:)'];
%     end
%     [Y I ] = max(W);
%     result = zeros(M,N);
%     for i = 1 : nPhases
%         v{i} = reshape( (I==i),[M,N]);
%         result = result+v{i}*c(i);
%     end
%     figure(5);
%      colormap(gray); hold on;
%      imshow(result,[]);
uold = u;
toc
end


% % Update c
% W = [];
% for i = 1 : nPhases
%     W = [W;u{i}(:)'];
% end
% [Y I ] = max(W);
% result = zeros(M,N);
% for i = 1 : nPhases
%     v= reshape( (I==i),[M,N]);
%     count = sum(sum(v(:)>0));
%     s = sum(v(:).*f(:));
%     c(i) = s/count;
% end

W = [];
for i = 1 : nPhases
    W = [W;u{i}(:)'];
end
[Y I ] = max(W);
result = zeros(M,N);
for i = 1 : nPhases
    u{i}= reshape( (I==i),[M,N]);
    result = result+u{i}*c(i);
end
figure(2);
imshow(result,[])

% figure;
% imagesc(u0,[0, 1]); axis off; axis equal;
% colormap(gray); hold on;
% contour(u{2}, [0,0], 'r');hold on;
end







